Finite-temperature properties of frustrated classical spins coupled to the lattice 
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We present extensive Monte Carlo simulations for a classical antiferromagnetic Heisenberg model 
with both nearest (Ji) and next-nearest (J2) exchange couplings on the square lattice coupled to 
the lattice degrees of freedom. The Ising-like phase transition, that appears for J2/J1 > 1/2 in 
the pure spin model, is strengthened by the spin-lattice coupling, and is accompanied by a lattice 
deformation from a tetragonal symmetry to an orthorhombic one. Evidences that the universality 
class of the transition does not change with the inclusion of the spin-lattice coupling are reported. 
Implications for Li2VOSi04, the prototype for a layered J1 — J2 model in the collinear regime, are 
also discussed. 

PACS numbers: 75.10.Hk, 75.40.Cx, 75.40.Mg 



I. INTRODUCTION 



The role of frustrating interactions in low-dimensional 
systems is a very important aspect in the modern the- 
ory of magnetism in solidsi and molecular clusters^ In 
particular, the presence of finite-temperature phase tran- 
sitions in two-dimensional systems with continuous spin- 
rotational symmetry, that are ruled out by a naive in- 
terpretation of the Mermin and Wagner theorem, has 
been clearly documented during the last years and con- 
tinues to attract much attention, because of the variety 
of phenomena that could be generated at low temper- 
ature. 3 ' 4 ' 5 ' 6 ' 7 ' 8,9 The main point is that the presence of 
frustrating interactions can induce non-trivial degrees of 
freedom, that may undergo a phase transition when the 
temperature is lowered. Among the strongly frustrated 
spin systems, probably one of the most important ex- 
amples is the antiferromagnetic Heisenberg model on the 
square lattice with both nearest (Ji) and next-nearest 
neighbor ( J2) couplings, for which it has recently become 
certain that a very interesting scenario shows up for large 
enough frustrating ratio J2I J\- Indeed, for J2/J1 > 1/2 
and classical spins, the two sublattices are completely 
decoupled at zero temperature, and each of them has 
an independent antiferromagnetic order. Therefore, the 
ground-state energy does not depend on the relative ori- 
entation between the magnetizations of the two sublat- 
tices and the ground state has an 0(3) x 0(3) symmetry. 
At low temperature, thermal fluctuations lift this huge 
degeneracy by an order by disorder mechanismAS. and two 
families of collinear states are entropically selected, with 
pitch vectors Q = (0,7r) and (ir, 0), respectively. Chan- 
dra, Coleman, and Larking argued that this fact reduces 
the symmetry to 0(3)xi?2, and that the Z2 symmetry 
is broken at low temperature, giving rise to an Ising-like 
phase transition at finite temperature. Recently^ it has 
been possible to verify this scenario by using extensive 
Monte Carlo simulations, and a rather accurate estimate 



of the critical temperature as a function of the ratio J2/ J\ 
has been obtained. 

The fact that the two states with Q = (0, it) and (ir, 0) 
break the rotational symmetry, having antiferromagnetic 
spin correlations in one spatial direction and ferromag- 
netic correlations in the other, suggests that, once the 
spin-lattice coupling is considered, the lattice could also 
experience a phase transition, ferromagnetic and antifer- 
romagnetic bonds acquiring different lengths. Indeed, on 
general grounds, the super-exchange couplings come from 
a virtual hopping of the electrons, which depends upon 
the actual distance of the ions. A first evidence that a 
similar distortion of the lattice appears when the spins 
are coupled to classical lattice distortions has been found 
in Rcf. I 111 where the quantum J1 — J2 model at zero tem- 
perature has been analyzed by Lanczos diagonalization 
on small clusters. In particular, by considering a 4 x 4 
lattice and spin-1/2 coupled to adiabatic phonons, it has 
been shown that there is a large region where the lattice 
distorts and ferromagnetic and antiferromagnetic bonds 
acquire different lengths. 

Interestingly, there is also clear evidence that such 
a scenario is realized in a real compound, Li2VOSi04, 
a vanadate which can be considered as a prototype of 
the J1—J2 model in the collinear regiontiSiH Indeed, al- 
though the value of J2/J1 is not exactly known 
all estimates indicate that J 2 > J\. Moreover, NMR and 
muon spin rotation magnetization provide a clear evi- 
dence for the presence of a phase transition to a collinear 
order at T c ~ 2.&K and a structural distortion at a 
nearby temperature. A simple and appealing explana- 
tion relies on the existence of the Ising-like transition 
and the concomitant lattice distortions. 

Since an unbiased finite-temperature analysis of the 
quantum model for large clusters is, at present, impos- 
sible, in this work we would like to address the simpler 
problem of classical spins coupled to lattice distortions. 
Although this could be viewed as a crude approximation 
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to the original quantum model, especially for low tem- 
perature and near the disordered region (expected near 
J2/J1 ~ 1/2), the classical case of spins coupled to lat- 
tice deformation represents a highly non-trivial problem 
and, certainly, gives the zeroth-order approximation for 
the true quantum case. 

The paper is organized as follows: In Sec. [HI we present 
the model and the method we used and in Sec. IIIII wc 
show the results and draw the conclusions. 



II. THE MODEL AND THE METHOD 

In this section we introduce the spin-lattice Hamilto- 
nian and briefly describe the method we used to treat the 
lattice degrees of freedom. The Hamiltonian reads: 

H — Ji(«%)Sj • Sj + ^2 J2(dij)Si ■ Sj + 



I<2 



d° 



"'j 



,(1) 



where S, are 0(3) spins on a periodic square lattice with 
N = L x L sites. and indicate the sum 

over nearest and next-nearest neighbors, respectively. 
The super-exchange couplings J\(dij) and J2{dij) depend 
upon the distance dij of the sites i and j, and K\ and K2 
are the elastic coupling constants. Finally, d\^ is the bare 
lattice distance. In transition metal compounds, general 
arguments lead to exchange integrals that vary like the 
inverse of the distance to a certain power 9, with 9 in the 
range 5— 15. 16 Therefore, for small displacements around 
the equilibrium positions, we can write: 



J(d ij ) = J[$L 




(2) 



In the Hamiltonian of Eq. QJ, we treat both the spins 
Sj and the lattice coordinates dij as dynamical variables, 
allowing them to change their configurations. With re- 
spect to the pure spin problem, this fact doubles the 
number of dynamical variables present in the problem 
and makes the Monte Carlo algorithm much heavier. In 
particular, we find better to use Monte Carlo algorithms 
with local and/or global updates, instead of more in- 
volved Monte Carlo methods based on the reconstruc- 
tion of the density of statesjiii^ that we used in the case 
where the lattice is kept fixed. 

In practice, we work on a torus with two independent 
lengths L x and L yi which play the role of additional de- 
grees of freedom sampled by Monte Carlo. Moreover, 
each site of the cluster can independently change its po- 
sition in the lattice. These facts allow us to consider 
structural distortions as well as the usual spin proper- 
ties. We found that, besides small displacements around 
their equilibrium positions, the sites always form either 
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FIG. 1: Zero-temperature phase diagram of the Hamiltonian 
of Eq. Q: the Neel state, with a tetragonal lattice, is stabi- 
lized for J2/J1 < 0.5 and below a critical value of the spin- 
lattice coupling 6 C , otherwise the ground state has collinear 
spins, with antiferromagnetic correlations along one spatial 
direction and ferromagnetic along the other, and the lattice 
has an orthorhombic structure. Lines are guides to the eye. 



a square cluster, with tetragonal symmetry and equal 
lattice spacing, l x and l y , in the two spatial directions 
(i.e., l x = l y ) or a rectangular cluster, with orthorhom- 
bic symmetry and different lattice spacings in the two 
directions (i.e., l x ^ l y ). In these cases the total lengths 
of the cluster in the two directions are L x = L x l x and 

Ly = L X ly . 

For the pure spin model, e.g., with fixed distances and 
J(dij) = J, in order to characterize the Ising-like phase 
transition, it is useful to construct, from the original spin 
variables Sj, an effective Ising variable (on the dual lat- 
tice)^ 



(Sj - Sfc) ■ (Sj- - SQ 

|(s< - s fc ) • fa - sol 



(3) 



where (i,j,k,l) are the corners with diagonal (i,k) and 
(j, I) of the plaquette centered at the site x of the dual lat- 
tice. In this way, the two collinear states with Q = (ir, 0) 
and Q = (0, n) can be distinguished by the value of the 
Ising variable, a x — ±1. As emphasized in Ref. 0, the 
normalization term does not affect the critical properties 
of the model and it is only introduced to have a normal- 
ized variable. 

Moreover, we can easily construct the Ising-like order 
parameter as M a — (^/N)^ x a x . Associated to the 
Ising magnetization M a , we can define the susceptibility 
related to the Ising variable % = {N/T)((M%) - (\M a \) 2 ). 
Finally, in order to study the finite-temperature phase di- 
agram, it is also useful to consider the specific heat per 
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siteCV, = ( 1 / LT 2 )({E 2 ) — {E) 2 ),E being the total energy 
of the system. 

Before considering the Monte Carlo results at finite 
temperature, it is useful to discuss the zero-temperature 
phase diagram of the spin-lattice Hamiltonian of Eq. , 
see Fig. ^ hi this case, the first-order phase transition 
present for 9 = 0, i.e., for the pure spin model, at J%l J\ = 
0.5 bends towards smaller values of the frustrating ratio 
J 2 /Ji when the coupling to the lattice is switched on. 
In particular, the huge 0(3) x 0(3) degeneracy is already 
broken at zero temperature, and the collinear states with 
pitch vector Q = (it, 0) and Q = (0, it) have a lower en- 
ergy. It is worth noting that this is a general property 
of the spin-lattice coupling, which is also present in the 
quantum case^i Therefore, on the basis of the analy- 
sis presented in Ref. Q, in the presence of the spin-lattice 
coupling, we expect that a finite-temperature phase tran- 
sition shows up also for a bare ratio J%/ J\ < 1/2, when- 
ever the coupling 9 is sufficiently large, i.e., 9 > 9 C . In the 
following, we will present strong evidences that there is 
a finite-temperature phase transition that is related to a 
change in the lattice symmetry, from a high-temperature 
disordered magnetic phase with a tetragonal lattice to 
a low-temperature Ising-ordered phase, with M a i= 0, 
and with an orthorhombic lattice. A little bit more sub- 
tle is the determination of the universality class of this 
transition. By using powerful Monte Carlo methods, we 
showed that, in the absence of spin-lattice coupling, the 
transition belongs to the Ising universality class^ as ex- 
pected from more general arguments^ The concomitant 
presence of spin and lattice degrees of freedom makes the 
Monte Carlo calculation much heavier than in the sim- 
pler case where only spins are considered, and we are 
limited to smaller clusters. Nevertheless, we can reach 
rather large lattices, that enable us to have convincing 
results for the thermodynamic limit. 



III. RESULTS AND DISCUSSION 

In this section, we present the Monte Carlo results for 
the Hamiltonian of Eq. and we give evidences that 
the Ising transition present in the pure spin model sur- 
vives when the spin-lattice coupling is taken into account, 
and, moreover, that it is accompanied by a structural de- 
formation. 

As discussed in the previous section, in the presence of 
both spin and lattice degrees of freedom, the Monte Carlo 
simulations become rather expensive and we are limited 
concerning the size of the clusters. In order to have a 
sufficiently high transition temperature T c , that allows a 
reliable estimate of the physical properties, most of the 
calculations have been performed for Jij J\ = 0.8 and 
three different values of the spin-lattice coupling, 9 = 5, 
10, and 15 for K\ = K2 = K. In the following we con- 
sider the case of KjJ\ = 100. Then, we also report 
calculations for other values of the ratio J2/J1, where a 
similar behavior is observed. Therefore, on the basis of 
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FIG. 2: Vmin(L) as a function of 1/L 2 for 9 = 5, 10, and 15. 
The horizontal line, located at 2/3, indicates the value for a 
continuous phase transition. 
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FIG. 3: Upper panel: Specific heat as a function of the 
temperature for 9 = 5 for three different clusters L = 20, 40, 
and 80. Lower panel: The same for the spin susceptibility. 



our numerical results, we argue that different choices for 
the frustrating ratio J%j J\ or for the elastic coupling K\ 
and K2 do not change qualitatively the physical behavior 
of the system. However, it should be emphasized that a 
precise determination of the critical exponents could be 
difficult and, in general, rather large clusters are needed. 
In the following, we present evidences that the transi- 
tion belongs to the Ising universality class, which has 
been clearly documented in the absence of spin-lattice 
coupling. 
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FIG. 4: The same as in Fig. for 9 = 10. 



First of all, in order to have some useful insight into 
the order of the phase transition, we consider Binder's 
fourth energy cumulant: 



V L (T) = 1 



(E 4 ) 
3(E 2 )' 



(4) 



Indeed, by considering the size scaling of V m in(L), the 
value of Vl {T) at the "critical" temperature of the finite 
cluster T C (L), it is possible to discriminate between a con- 
tinuous transition, for which V* = limi-,^, V m i n (L) = 
2/3, and a first-order transition, for which V* < 2/3. 19 
The size dependence of Vl(T c ) is shown in Fig. [21 for 
three values of 9. Interestingly, in all the cases consid- 
ered V m i rl (L) — > 2/3, giving a clear evidence that the 
transition is continuous and not first order. 

In Figs. and 0] we show the specific heat and the 
spin susceptibility x f° r ^ — 5 and 10, respectively. In 
both quantities, the presence of a huge peak, that grows 
with the size of the cluster, marks the presence of a phase 
transition. It is worth noting that, while the maximum 
value of the spin susceptibility is not affected by the spin- 
lattice coupling, the maximum value of the specific heat 
strongly depends upon 9. Moreover, the low-temperature 
value of the specific heat is renormalized by the presence 
of the lattice displacements: the bare value C v = 1, that 
holds when the lattice is frozen in its equilibrium config- 
uration, is changed into C v = 2 whenever the sites can 
have small displacements around their equilibrium posi- 
tions. 

The actual determination of the universality class of 
the transition can be assessed by considering the size scal- 
ing of the physical quantities. For a second-order phase 
transition, the maximum value of the susceptibility and 
of the specific heat follow well defined scaling relations, 
i.e., Xmax(L) ~ LPl v and C max (L) ~ L a (a logarithmic 
behavior, i.e., a = 0, for the Ising universality class). 



On the other hand, for a discontinuous phase transition, 
we have Xmax(L) ~ L 2 and C max (L) ~ L 2 . As clearly 
shown for 9 = 0j£ when the spins are completely decou- 
pled from the underlying lattice, the phase transition is 
second order and falls into the Ising universality class, 
i.e., a = 0, v = 1 and 7 = 7/4. In this case, an ac- 
curate determination of a is quite hard and rather large 
clusters are needed (i.e., L ~ 200). In the case of a 
finite 9, we cannot afford such large clusters, but, never- 
theless, we can obtain rather convincing evidences in fa- 
vor of the Ising universality class by considering the case 
of a large spin-lattice coupling 9. Indeed, in Fig. [SJ we 
show the size scaling of the maximum of the specific heat 
and of the spin susceptibility for 9—15. As mentioned 
before, the value of the peak of the spin susceptibility 
does not depend upon 9, and, therefore, the results are 
also valid for other values of the spin-lattice couplings. 
A three-parameter fit of Xmax(L) = a + b x L 1 ^ gives 
7/1/ = 1.7±0.1, which is in reasonable agreement with the 
value 1.75 of the Ising transition. Moreover, the clearest 
evidence that the transition belongs to the Ising univer- 
sality class comes from the specific heat. Indeed, the 
best fit of the maximum of the specific heat is given by 
Cmax(L) — ciq + a\ log(L) + a2/L^ while a power law 
turns out to be completely inadequate. This fact com- 
pletely rules out the possibility to have a first-order phase 
transition, for which a power law with an exponent a = 2 
is expected. For smaller values of 9 it is much harder to 
get an accurate fit of C max (L\ and much larger sizes 
should be considered. Nevertheless, having in mind the 
results for 9 = 0, we expect that the phase transition also 
falls into the Ising universality class for small spin-lattice 
couplingiSi 

By using a similar analysis for different J2/J1, it is 
possible to determine the behavior of the critical temper- 
ature T c as a function of the frustrating ratio Ji / J\ for 
different values of 9. The results are reported in Fig. 
In particular, in order to determine T c , we used two al- 
ternative methods that give consistent results: On one 
hand, we take the temperatures for which the suscepti- 
bility has its maximum and we make use of the scaling 
T C (L) ~ T c + ax L^ 1 (which assumes v = 1); on the other 
hand, we used Binder's fourth cumulant of the Ising pa- 
rameter. Interestingly, as already anticipated, for a suf- 
ficiently large value of the spin-lattice coupling there is a 
phase transition to a collinear phase even if J2/J1 < 1/2 
as soon as the spin-lattice coupling is switched on, the 
critical ratio J2/J1 for which a transition occurs for a 
given 9 being in good agreement with the static analysis 

(Fig. nj. 

Now, we turn to the discussion of the lattice proper- 
ties. Because the two families of states with pitch vec- 
tors Q = (0, 7r) and (tt,0), that are entropically selected, 
have different spin correlations in the two spatial direc- 
tions, whenever the sites of the lattices are coupled to 
the spins, a structural transition is also likely to show 
up. Indeed, while at high temperature the lattice has a 
tetragonal symmetry, with the same lattice spacing in the 
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FIG. 5: Upper panel: Size scaling of the maximum of the 
specific heat as a function of the size of the cluster for 9 = 15. 
Lower panel: The same for the maximum of the spin suscep- 
tibility. In both cases, the dashed line is a three-parameter 
fit (see text). 



FIG. 7: Upper panel: lattice parameters l x and l y as a func- 



tion of the temperature for 9 = 5 for L — 20 and 80. 
panel: The same for the Ising magnetization M a . 
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FIG. 6: The critical temperature T c as a function of the frus- 
trating ratio J2/J1 for three different values of the spin-lattice 
couplings, = 5, 10, and 15. The case for the pure spin model, 
corresponding to 6 = 0, is also reported for comparison. 



x and y directions, at the transition temperature, the lat- 
tice undergoes a structural phase transition towards an 
orthorhombic symmetry, with different lattice spacings 
l x and l y , see Figs. [7] and 00 This is exactly equivalent to 
what happens at zero temperature in the quantum case, 
where the collinear phase is also accompanied by a lattice 
distortion. Notice that, since we are considering classical 
spins, the only relevant distortion is the orthorhombic 



FIG. 8: The same as in Fig. for 6 = 10. 



one, and, indeed, we do not find any evidence towards 
other displacements of the underlying lattice, in contrast 
to the spin 1/2 caseiii 

One important outcome is that the structural phase 
transition has a huge effect just below the critical tem- 
perature. Indeed, the lattice displacements strongly sta- 
bilize one of the two families of states, and almost freeze 
the spin fluctuations. This can be clearly seen from the 
Ising-like magnetization M a , whose value is practically 
fixed to its saturation value from zero temperature up to 
T c , where it goes to zero (or to a very small value for 
finite sizes), see Figs. [7| and 00 
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The implications for quasi-two-dimensional models 
with inter-plane magnetic couplings will presumably de- 
pend on the value of this coupling. If it is sufficiently 
weak, we expect the system to undergo two phase transi- 
tions: First, coming from high temperature, an Ising like 
transition involving a lattice distortion and a breaking 
of the rotational symmetry of spin fluctuations, then a 
three-dimensional magnetic ordering. However, since the 
correlation length increases exponentially fast at low tem- 
perature in two-dimensional antifcrromagnets, the crit- 
ical temperature for three-dimensional ordering is ex- 
pected to grow very fast with the inter-plane coupling, 
reaching quite rapidly values of the order of the intra- 
plane couplings. In that situation, the two transitions 
might very well merge into a single phase transition, with 
a critical behavior that could be very different from the 
standard one for magnetic ordering of three-dimensional 
Heisenberg antifcrromagnets though.— Given the lattice 
sizes available when phonons arc included, this is an issue 
we could not address however. 

From an experimental point of view, the first proto- 
type of the spin- 1/2 J1 — J2 model has been recently syn- 
thetized^ It is a layered vanadium oxide Li2VOSi04 in 
which VO5 pyramids are arranged in such a way that sec- 
ond vanadium neighbors are in the same plane while first 
neighbors are not, so that J2 need not be smaller than 
Ji . Although there is no general agreement on the precise 
value of the ratio J2/ J\ , both experimental and theoreti- 
cal estimates lead to a ratio significantly larger than 1/2. 
Therefore, the system is expected to develop collinear 
order, which has been confirmed by NMR results at the 
Li site. Indeed, the NMR spectrum at this site clearly 
shows that, below T c ~ 2.8K, the central peak splits into 



three different peaks that correspond to different Li sites 
where the local field is either zero or nonzero (parallel or 
antiparallel to the external field). This behavior has been 
associated to a magnetic order where the spins lie along 
the x direction with the magnetic vector Q = (0, 7r)^£ 
Moreover, the 29 Si NMR spectrum is very anomalous at 
low temperature^ and given the very symmetric posi- 
tion of Si in the lattice, this cannot be attributed to the 
development of collinear order but must be related to a 
structural phase transition. Given the rather low char- 
acteristic temperatures in that system (T c ~ 2.&K while 
the transfer of weight of the Si line is complete at 2 K), a 
precise structural determination of the low temperature 
phase has not been possible yet. Interestingly enough, 
the temperature at which the system starts to develop a 
structural distortion seems to be slightly larger than T Cl 
which would agree with the general expectation for small 
inter-plane coupling. Indeed, experimentally, it is not 
possible to exclude that distorted regions start to gener- 
ate above 2.8K and grow, by decreasing the temperature, 
by a nucleation procedure^ Whether a more direct de- 
termination of the structural transition will confirm this 
point remains to be seen. 
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